/*=auto=========================================================================

  Portions (c) Copyright 2005 Brigham and Women's Hospital (BWH) All Rights Reserved.

  See Doc/copyright/copyright.txt
  or http://www.slicer.org/copyright/copyright.txt for details.

  Program:   3D Slicer
  Module:    $RCSfile: vtkFSSurfaceReader.cxx,v $
  Date:      $Date: 2006/08/04 18:48:03 $
  Version:   $Revision: 1.14 $

=========================================================================auto=*/
#include "vtkFSSurfaceReader.h"
#include "vtkObjectFactory.h"
#include "vtkByteSwap.h"
#include "vtkFloatArray.h"
#include "vtkFSIO.h"
#include "vtkPolyData.h"
#include "vtkCellArray.h"
#if (VTK_MAJOR_VERSION >= 5)
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkInformation.h"
#endif

// -------------------------------------------------------------------------
vtkFSSurfaceReader * vtkFSSurfaceReader::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkFSSurfaceReader");

  if( ret )
    {
    // return (vtkFSSurfaceReader*)ret;
    return dynamic_cast<vtkFSSurfaceReader *> ret;
    }
  // If the factory was unable to create the object, then create it here.
  return new vtkFSSurfaceReader;
}

// -------------------------------------------------------------------------
vtkFSSurfaceReader::vtkFSSurfaceReader()
{
  vtkPolyData *output = vtkPolyData::New();

  this->SetOutput(output);

  // Releasing data for pipeline parallism.
  // Filters will know it is empty.
  output->ReleaseData();
  output->Delete();
  this->ExecutePiece = this->ExecuteNumberOfPieces = 0;
  this->ExecuteGhostLevel = 0;
}

// -------------------------------------------------------------------------
vtkFSSurfaceReader::~vtkFSSurfaceReader()
{
}

// ----------------------------------------------------------------------------
vtkPolyData * vtkFSSurfaceReader::GetOutput()
{
  int numberOfOutputs;

#if (VTK_MAJOR_VERSION >= 5)
  numberOfOutputs = this->GetNumberOfOutputPorts();
#else
  numberOfOutputs = this->NumberOfOutputs;
#endif

  if( numberOfOutputs < 1 )
    {
    return NULL;
    }
#if (VTK_MAJOR_VERSION >= 5)
  return dynamic_cast<vtkPolyData *>(this->GetOutput(0) );
#else
  return dynamic_cast<vtkPolyData *>(this->Outputs[0]);
#endif
}

// ----------------------------------------------------------------------------
vtkPolyData * vtkFSSurfaceReader::GetOutput(int idx)
{
#if (VTK_MAJOR_VERSION >= 5)
  return vtkPolyData::SafeDownCast( this->GetOutputDataObject(idx) );
#else
  return vtkPolyData::SafeDownCast( this->vtkSource::GetOutput(idx) );
#endif
}

// ----------------------------------------------------------------------------
void vtkFSSurfaceReader::SetOutput(vtkPolyData *output)
{
#if (VTK_MAJOR_VERSION >= 5)
  this->GetExecutive()->SetOutputData(0, output);
#else
  this->vtkSource::SetNthOutput(0, output);
#endif
}

// ----------------------------------------------------------------------------
#if !(VTK_MAJOR_VERSION >= 5)
void vtkFSSurfaceReader::ComputeInputUpdateExtents(vtkDataObject *data)
{
  int          piece, numPieces, ghostLevel;
  vtkPolyData *output = (vtkPolyData *)data;
  int          idx;

  output->GetUpdateExtent(piece, numPieces, ghostLevel);

  // make sure piece is valid
  if( piece < 0 || piece >= numPieces )
    {
    return;
    }

  if( ghostLevel < 0 )
    {
    return;
    }
  // just copy the Update extent as default behavior.
  for( idx = 0; idx < this->NumberOfInputs; ++idx )
    {
    if( this->Inputs[idx] )
      {
      this->Inputs[idx]->SetUpdateExtent(piece, numPieces, ghostLevel);
      }
    }

  // Save the piece so execute can use this information.
  this->ExecutePiece = piece;
  this->ExecuteNumberOfPieces = numPieces;

  this->ExecuteGhostLevel = ghostLevel;
}

#endif

int vtkFSSurfaceReader::RequestData(
  vtkInformation *,
  vtkInformationVector * *,
  vtkInformationVector *outputVector)
{
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
//  vtkPolyData *output = this->GetOutput();
  vtkPolyData *output = vtkPolyData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT() ) );
  FILE*         surfaceFile;
  int           magicNumber;
  char          line[256];
  int           numVertices = 0;
  int           numFaces = 0;
  int           vIndex, fIndex;
  int           numVerticesPerFace = 0;
  int           tmpX, tmpY, tmpZ;
  int           tmpfIndex = 0;
  float         locations[3];
  int           fvIndex;
  int           faceIncrement = 1;
  int           faceMultiplier = 1;
  vtkIdType     faceIndices[4];
  vtkPoints *   outputVertices;
  vtkCellArray *outputFaces;

#if FS_CALC_NORMALS
  vtkFloatArray *outputNormals;
  FSVertex*      vertices;
  FSFace*        faces;
  FSVertex*      v;
  FSFace*        f;
  FSVertex*      fv0;
  FSVertex*      fv1;
  FSVertex*      fv2;
  int            fvIndex1, fvIndex2, fvIndex0;
  float          faceVector0[3], faceVector1[3];
  float          faceNormal[3];
  float          length;
#endif
  int totalSteps = 1;
  int thisStep = 0;

  vtkDebugMacro(<< "RequestData: Reading vtk polygonal data...");

  // Try to open the file.
  surfaceFile = fopen(this->FileName, "rb");
  if( !surfaceFile )
    {
    vtkErrorMacro(<< "Could not open file " << this->FileName);
    return 1;
    }

  // Get the three byte magic number. We support two file types.
  vtkFSIO::ReadInt3(surfaceFile, magicNumber);
  if( magicNumber != vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER &&
      magicNumber != vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER &&
      magicNumber != vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER )
    {
    vtkErrorMacro(
      << "vtkFSSurfaceReader.cxx Execute: Wrong file type when loading " << this->FileName << "\n magic number = " << magicNumber << ". Supported ar " << vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER << ", " << vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER << ", and "
      << vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER );
    return 1;
    }

#if FS_DEBUG

  switch( magicNumber )
    {
    case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
      {
      cerr << "Reading old quad file" << endl;
      }
      break;
    case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
      {
      cerr << "Reading new quad file" << endl;
      }
      break;
    case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
      {
      cerr << "Reading triangle file" << endl;
      }
      break;
    }
#endif

  // Triangle file has some kind of header string at the
  // beginning. Skip it.
  if( vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER == magicNumber )
    {
    char *skipchars = fgets(line, 200, surfaceFile);
    int   skip = fscanf(surfaceFile, "\n");
    if( skipchars == NULL || skip > 0 )
      {
      // trying to avoid unused var warnings while checking return values
      }
    }

  // Triangle files use normal ints to store their number of vertices
  // and faces, while quad files use three byte ints.
  size_t retval;

  switch( magicNumber )
    {
    case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
    case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
      {
      vtkFSIO::ReadInt3(surfaceFile, numVertices);
      vtkFSIO::ReadInt3(surfaceFile, numFaces);
      }
      break;
    case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
      {
      retval = fread(&numVertices, sizeof(int), 1, surfaceFile);
      if( retval == 1 )
        {
        vtkByteSwap::Swap4BE(&numVertices);
        }
      else
        {
        vtkErrorMacro("Error reading number of vertices");
        }
      retval = fread(&numFaces, sizeof(int), 1, surfaceFile);
      if( retval == 1 )
        {
        vtkByteSwap::Swap4BE(&numFaces);
        }
      else
        {
        vtkErrorMacro("Error reading number of faces");
        }
      }
      break;
    }

  // In quad files, we want to skip every other face but count twice
  // as many of them. This has to do with the way they are stored;
  // here we just generate quads where as in the old code they
  // generated tries from the quads. (Trust me.) In tri files, we use
  // every face.
  switch( magicNumber )
    {
    case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
    case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
      {
      faceIncrement = 2;
      faceMultiplier = 2;
      }
      break;
    case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
      {
      faceIncrement = 1;
      faceMultiplier = 1;
      }
      break;
    }

#if FS_DEBUG
  cerr << numVertices << " vertices, " << numFaces * faceMultiplier << " faces" << endl;
#endif

  // If quad files, there are four vertices per face, in tri files,
  // there are three.
  switch( magicNumber )
    {
    case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
    case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
      {
      numVerticesPerFace = vtkFSSurfaceReader::FS_NUM_VERTS_IN_QUAD_FACE;
      }
      break;
    case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
      {
      numVerticesPerFace = vtkFSSurfaceReader::FS_NUM_VERTS_IN_TRI_FACE;
      }
      break;
    }

  // Allocate our VTK arrays.
  outputVertices = vtkPoints::New();
  outputVertices->Allocate(numVertices);
  outputFaces = vtkCellArray::New();
  outputFaces->Allocate(outputFaces->EstimateSize(numFaces,
                                                  numVerticesPerFace) );
#if FS_CALC_NORMALS
  outputNormals = vtkFloatArray::New();
  outputNormals->Allocate(numVertices);
  outputNormals->SetNumberOfComponents(3);
  outputNormals->SetName("Normals");

  // Allocate our vertex and face connectivity arrays for calculating
  // the normals. If we can't, no big deal.
  vertices = (FSVertex *) calloc(numVertices, sizeof(FSVertex) );
  faces = (FSFace *) calloc(numFaces, sizeof(FSFace) );
  if( NULL == vertices || NULL == faces )
    {
    vtkErrorMacro(<< "Couldn't allocate vertex or face connectivity "
                  "structures, unable to calc normals");
    }
#endif

  totalSteps = numVertices + (numFaces * faceMultiplier / faceIncrement * numVerticesPerFace);
#if FS_CALC_NORMALS
  if( NULL != vertices && NULL != faces )
    {
    for( vIndex = 0; vIndex < numVertices; vIndex++ )
      {
      fv1 = &vertices[vIndex];
      for( fvIndex = 0; fvIndex < fv1->numFaces; fvIndex++ )
        {
        totalSteps++;
        }
      }
    }
#endif
  vtkDebugMacro(<< "Got total steps = " << totalSteps);
  // For each vertex...
  for( vIndex = 0; vIndex < numVertices; vIndex++ )
    {
    thisStep++;

    // Depending on the file type, read in three two bytes ints and
    // convert them from meters to millimeters or read in three floats
    // in millimeters. Insert them into the vertices array.  The old
    // quad format uses the ints and the new quad and triangle formats
    // use floats.
    switch( magicNumber )
      {
      case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
        {
        vtkFSIO::ReadInt2(surfaceFile, tmpX);
        vtkFSIO::ReadInt2(surfaceFile, tmpY);
        vtkFSIO::ReadInt2(surfaceFile, tmpZ);
        locations[0] = (float)tmpX / 100.0;
        locations[1] = (float)tmpY / 100.0;
        locations[2] = (float)tmpZ / 100.0;
        }
        break;
      case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
      case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
        {
        vtkFSIO::ReadFloat(surfaceFile, locations[0]);
        vtkFSIO::ReadFloat(surfaceFile, locations[1]);
        vtkFSIO::ReadFloat(surfaceFile, locations[2]);
        }
        break;
      }

    outputVertices->InsertNextPoint(locations);

#if FS_CALC_NORMALS
    // If we have connectivity info, fill out the location and
    // initialize the normals to 0 for this vertex.
    if( NULL != vertices && NULL != faces )
      {
      v = &vertices[vIndex];
      v->x = locations[0];
      v->y = locations[1];
      v->z = locations[2];
      v->nx = 0;
      v->ny = 0;
      v->nz = 0;
      v->numFaces = 0;
      }

#endif
    if( (thisStep % 1000) == 0 )
      {
      this->UpdateProgress(1.0 * thisStep / totalSteps);
      }
    }
  // For each face...
  for( fIndex = 0;
       fIndex < numFaces * faceMultiplier;
       fIndex += faceIncrement )
    {
    // For each vertex in the face...
    for( fvIndex = 0; fvIndex < numVerticesPerFace; fvIndex++ )
      {
      thisStep++;

      // Read in a vertex index. Triangle format gets a normal int,
      // quad formats get three byte ints.
      switch( magicNumber )
        {
        case vtkFSSurfaceReader::FS_QUAD_FILE_MAGIC_NUMBER:
        case vtkFSSurfaceReader::FS_NEW_QUAD_FILE_MAGIC_NUMBER:
          {
          vtkFSIO::ReadInt3(surfaceFile, tmpfIndex);
          }
          break;
        case vtkFSSurfaceReader::FS_TRIANGLE_FILE_MAGIC_NUMBER:
          {
          retval = fread(&tmpfIndex, sizeof(int), 1, surfaceFile);
          if( retval == 1 )
            {
            vtkByteSwap::Swap4BE(&tmpfIndex);
            }
          else
            {
            vtkErrorMacro("Error reading integer at index " << fvIndex);
            }
          }
          break;
        }

      faceIndices[fvIndex] = tmpfIndex;

#if FS_CALC_NORMALS
      // Fill out connectivity info for this face. Get the vertex from
      // the vertex index and add this face index to its list of
      // faces. Then add this vertex index to the list of indicies in
      // the face.
      if( NULL != vertices && NULL != faces )
        {
        v = &vertices[tmpfIndex];
        v->faces[v->numFaces] = fIndex;
        v->indicesInFace[v->numFaces] = fvIndex;
        v->numFaces++;

        f = &faces[fIndex];
        f->vertices[fvIndex] = tmpfIndex;
        }
#endif
      }

    // Add the face to the list.
    outputFaces->InsertNextCell(numVerticesPerFace, faceIndices);
    if( (thisStep % 1000) == 0 )
      {
      this->UpdateProgress(1.0 * thisStep / totalSteps);
      }
    }

  // Close the surface file.
  fclose(surfaceFile);

#if FS_DEBUG
  cerr << "Done reading surface." << endl;
#endif

#if FS_CALC_NORMALS

  // If we allocated the space for our vertex and face connectivity,
  // calculate the normals.
  if( NULL != vertices && NULL != faces )
    {
    // For each vertex...
    for( vIndex = 0; vIndex < numVertices; vIndex++ )
      {
      // For each face it is a part of...
      fv1 = &vertices[vIndex];
      for( fIndex = 0; fIndex < fv1->numFaces; fIndex++ )
        {
        thisStep++;

        // Get this face. fv1->indicesInFace tells us which index,
        // 0 - numVerticesPerFace-1, it is in the face. Get the two
        // indicies surrounding it so we can get the vertices
        // surrounding it. Now we have fv0, fv1, and fv2 which are all
        // adjacent in the face, with fv1 being the orignal vertex.
        f = &faces[fv1->faces[fIndex]];
        fvIndex1 = fv1->indicesInFace[fIndex];
        if( 0 == fvIndex1 )
          {
          fvIndex0 = numVerticesPerFace - 1;
          }
        else
          {
          fvIndex0 = fvIndex1 - 1;
          }
        if( numVerticesPerFace == fvIndex1 )
          {
          fvIndex2 = 0;
          }
        else
          {
          fvIndex2 = fvIndex1 + 1;
          }

        fv0 = &vertices[f->vertices[fvIndex0]];
        fv2 = &vertices[f->vertices[fvIndex2]];

        // Get two vectors from these points and normalize them. Then
        // just cross them to get the perpendicular vector. This is
        // the normal for the face. Add this to the normal for the
        // vertex.
        faceVector0[0] = fv1->x - fv0->x;
        faceVector0[1] = fv1->y - fv0->y;
        faceVector0[2] = fv1->z - fv0->z;
        faceVector1[0] = fv2->x - fv1->x;
        faceVector1[1] = fv2->y - fv1->y;
        faceVector1[2] = fv2->z - fv1->z;
        length = sqrt( faceVector0[0] * faceVector0[0]
                       + faceVector0[1] * faceVector0[1]
                       + faceVector0[2] * faceVector0[2] );
        if( length > 0 )
          {
          faceVector0[0] /= length;
          faceVector0[1] /= length;
          faceVector0[2] /= length;
          }
        length = sqrt( faceVector1[0] * faceVector1[0]
                       + faceVector1[1] * faceVector1[1]
                       + faceVector1[2] * faceVector1[2] );
        if( length > 0 )
          {
          faceVector1[0] /= length;
          faceVector1[1] /= length;
          faceVector1[2] /= length;
          }

        // get the cross product
        faceNormal[0] = -faceVector1[1] * faceVector0[2]
          + faceVector0[1] * faceVector1[2];
        faceNormal[1] = faceVector1[0] * faceVector0[2]
          - faceVector0[0] * faceVector1[2];
        faceNormal[2] = -faceVector1[0] * faceVector0[1]
          + faceVector0[0] * faceVector1[1];

        // add it to the normal vector at this vertex.
        fv1->nx += faceNormal[0];
        fv1->ny += faceNormal[1];
        fv1->nz += faceNormal[2];
        }

      // When all the faces for this vertex have been processed, the
      // normal at the vertex is the sum of all the normals for the
      // adjacent faces. Normalize it and we're done.
      faceNormal[0] = fv1->nx;
      faceNormal[1] = fv1->ny;
      faceNormal[2] = fv1->nz;
      length = sqrt( faceNormal[0] * faceNormal[0]
                     + faceNormal[1] * faceNormal[1]
                     + faceNormal[2] * faceNormal[2] );
      if( length > 0 )
        {
        faceNormal[0] /= length;
        faceNormal[1] /= length;
        faceNormal[2] /= length;
        }

      if( (thisStep % 1000) == 0 )
        {
        this->UpdateProgress(1.0 * thisStep / totalSteps);
        }

      // Add the final normal to the array.
      outputNormals->InsertNextTuple(faceNormal);
      }

    free(vertices);
    free(faces);
    }
#endif

  // Set all the arrays in the output.
  output->SetPoints(outputVertices);
  outputVertices->Delete();

  this->SetProgressText("");
  this->UpdateProgress(0.0);

#if FS_CALC_NORMALS
  output->GetPointData()->SetNormals(outputNormals);
  outputNormals->Delete();
#endif

  outputFaces->Squeeze();
  output->SetPolys(outputFaces);
  outputFaces->Delete();

  return 1;
}

// ----------------------------------------------------------------------------
#if (VTK_MAJOR_VERSION >= 5)
int vtkFSSurfaceReader::FillOutputPortInformation(int,
                                                  vtkInformation* info)
{
  info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
  return 1;
}

#endif

// ----------------------------------------------------------------------------
void vtkFSSurfaceReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
}
